Symbolic software for separation of variables in the 
Hamilton-Jacobi equation for the L-systems. 



Yu.A. Grigoryev and A.V. Tsiganov 



Department of Computational Physics, 
St. Petersburg State University, St. Petersburg, Russia 
e-mail: tsiganov@mph.phys.spbu.ru 

Abstract 

We discuss computer implementation of the known algorithm of finding separation 
coordinates for the special class of orthogonal separable systems called L-systems or Benenti 
systems. 

1 Introduction 

Let Q be an-dimensional Riemannian manifold with generic local coordinates q = (q , q 2 , . . . , q n ) 
and positive-definite metric tensor G. 

On the cotangent bundle T* Q of Q with canonical coordinates (p, q) we consider dynamical 
system with the natural Hamilton function 

n 

H = T{p,q) + V(q) = J2 ^(l)PiPo + V(q). (1.1) 

i,3=l 

Here g u (g) are components of the metric tensor G and V{q) is the potential energy, a smooth 
function on Q canonically lifted to a function on T*Q. 

One of the most effective ways of solving the corresponding equations of motion is by 
separation of variables in the Hamilton-Jacobi equation 

H{ Pl q) = E. (1.2) 

A coordinate system Q = (Q 1 , . . . , Q n ) is called separable if the Hamilton-Jacobi equation 
admits a complete solution of the form 



S{Q,a) = ^<S;(Q\a), det 



d 2 S 
dCfda? 



?0. 



i=l 

Here a = (or, . . . , a") is a set of separation constants. The corresponding Jacobi equations 

dSi{Q\a) 



are called the separated equations. In a similar manner, we will call a natural Hamiltonian or 
potential separable if such a separable coordinate system exists. 

For a given Hamilton function H(p, q) the problem of finding canonical transformation from 
initial variables (p, q) to the separation coordinates (P, Q) is very non-trivial. The problem was 
originally formulated by Jacobi when he invented elliptic coordinates and successfully applied 
them to solve several important mechanical problems, such as the problem of geodesic motion 
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on an ellipsoid and the Eulcr problem of planar motion in a force field of two attracting centers 
and the problem of the motion of three particles, which interact due to forces depending on 
their relative distances pQ. 

Up to now finding the separation coordinates for a given integrable system remains rather 
the magic art than a constructive theory. However, for the special class of the natural Hamil- 
tonians we have a complete and algorithmic solution of this problem 03 E] . The algorithm 
is straightforward enough to be implemented on a computer and thus turned into a practical 
tool. 

In this note we present an implementation of this algorithm by means of the computer 
algebra system Maple 9.5. The corresponding Maple file with the code may be found in 

http : //www.maplesof t . com/applications/app_center_view. aspx?AID=1686 

2 Algorithm of the point separation of variables 

According to |S] we restrict ourselves to the search of the point canonical transformations 
Q = /(<?) and P = g(q,p) only. 

In this case from Stackel (1893), Levi-Civita (1904), Eisenhart (1934), Kalnins & Miller 
(1980) and Benenti (1993) we have 

Theorem 1 The Hamilton- Jacobi equation fl.ty) is separable in orthogonal coordinates if and 
only if there exists symmetric Killing 2-tensor K with simple eigenvalues and normal eigenvec- 
tors such that 

d(KdV) = 0, (2.3) 

where d denotes the exterior derivative. 

Such Killing tensor K is called characteristic, its existence is completely defined by the kinetic 
energy T. It means that the separation of the geodesic equation is a necessary condition for 
the separation of equation T + V = E. The equation d{KdV) — is an integrability condition 
for the existence of the potential V(q) that may be added to T. 

In 1992 Benenti has shown a simple recurrence procedure to construct a special family of 
Killing tensors K obeying the assumptions of this theorem. He considered a special class of Rie- 
mannian manifolds Q endowed with the L-tensor, whose functionally independent eigenvalues 
are identified with the separation variables. 

Following Benenti let us call L-tensor a conformal Killing tensor L with vanishing torsion 
and pointwise simple eigenvalues Qi. Under these conditions the tensors 

m 

K m = ^2 cr TO _ fe L fc , or K m = a m G - K m _iL, m = 0, . . . , n - 1, (2.4) 

fc=0 

where L is (2,0) tensor field, are the Killing tensors with simple eigenvalues and normal eigen- 
vectors. Here functions the elementary symmetric polynomials of degree m on the 

n 

eigenvalues of L, such that det(AI — L) = er m A ,l ~ m . 

m=0 

It can be shown |S] that equation d(KdV) = implies d(K m dV) = for all the commuting 
Killing tensors. This allows to define the new potentials V m according to 

dV m = K m dV, m = 1, . . . , n - 1 (2.5) 

and to introduce integrals of motion 

n 

H m = J2 K 'i>l -I>.l>, + V m (q), m = 1, . . . , n - 1. (2.6) 
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These functions form a family of commuting integrals of the motion for the Hamiltonian H 
According to |U] integrals of motion H m l|2.6|) are solutions of the recursion relations 

dH m+1 = N*dH m + <r m+1 dH, m = 1, . . . , n - 1, H n = . (2.7) 
Here N is the recursion operator, which is the complete lifting of (1,1) tensor field L to T*Q 

dq k k g qi T y Qqk g qi J Qp , igpi - 

Recall, since a metric tensor is present, the boldface object K or L can be represented in 
components as a tensor of type (2, 0), (1, 1) and (0, 2), respectively. 

So, in order to get separation variables and integrals of motion for a given integrable system 
we have to find tensor L only. As noticed in [8], tensor L is an L-tensor with respect to the 
usual Riemannian metric G iff 

d(Cx T - Tdai) = 0, (2.8) 

where C is the Lie derivative along the geodesic vector field Xt, o\ — trL is first symmetric 
polynomial and 

n 

II ^ /.>,</</' (2.9) 

is the L-dcformation of the standard Liouville 1-form 8q — ^Pjdq* for any set of fibered 
coordinates (p, q) . 

The equation d(KdV) = with the Killing tensor Ki lf2~H> may be rewritten in the similar 

form 

d{C Xv e-Vda x ) =0, (2.10) 

see for instance [TU] . 

Below we will use the following well-known expression for the Lie derivative C along the 
vector field X 

C x = ±x d + d± x ■ 

Here ±x is a hook operator and d is an exterior derivative. It allows us to decrease the amount 
of intermediate calculations as long as d 2 = and we have 

dCx = d±x d + d 2 ±x = dix d. 

In these notations equations l|2.8|l and H2.10JI read as 

d(±x T d9 -Tda x ) = 0, (2.11) 
d(i Xv d0 -Vda-i) = 0. (2.12) 

We will call L-system or Benenti system any separable orthogonal system whose Killing 
tensor K in H2.3(l is generated by an L-tensor according to (|2.4|l . To construct the separation 
coordinates Q for the L-system we have to solve equations (I2.11f) and l|2.12|l with respect to 
functions L*(q) and to find the eigenvalues of the tensor L. 

In the next section we present a computer program for search of the L-tensors and the 
associated separation coordinates for L-systems. 

Remark: The 1-form 9 may be used to construct the second Poisson structure on T*Q [9"lll0|. 
In this case the recursion operator N will be a complete lifting of the L-tensor L from the 
configuration space to the whole phase space. Therefore, the coordinate separation of variables 
can be treated as a particular case of the bi-hamiltonian theory of separation of variables. 
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3 The program for the search of separation variables 



In this section we present an implementation of the discussed algorithm made in the symbolic 
computational system Maple v. 9. 5. 

We will use the standard Maple package liesymm, which was designed for construction of 
the differential forms corresponding to partial differential equations. In fact we are doing the 
inverse procedure, i.e. starting with the differential forms (|2.11I2.12|) we have to get a system 
of partial differential equations. 

Let us start with the following command 

> with (liesymm) : 

which makes the short form names of the functions of a Maple package available at the inter- 
active level. 

At the first step we need to determine dimension of a given configuration space Q and to 
suppose that the phase space T*Q is equipped with some canonical coordinates q — (q , . . . , q n ) 
and p = (pi, . . . ,pn): 

> n := 2; 

n := 2 

> q: =seq(q I I i , i=l . .n) : p : =seq(p I I i , i=l . .n) : var:=q,p: 

> setup (var) ; 

[ql, q2, pi, p2] 

Now we will look for a tensor field L with vanishing torsion and functionally independent 
eigenvalues. For a given integrable system tensor L has to satisfy equations (|2.11|) and H2.12fl 
which include first symmetric polynomial o\ on the eigenvalues of L 

> Sigma: =add(L [i , i] (q) , i=l . .n) ; 

cr:=L M (gl, g2)+L 2 ,2(gl, «2) 
and L-deformation 9 (|2.9|) of the standard Liouville form 

> theta:=add(add(L[i, j] (q)*p| |i*d(q| I j) ,i=l. .n) ,j=l. .n) ; 

9;=L 1:1 (ql, q2)pld(ql)+L 2>1 (ql, q2)p2d(ql) 
+Li >2 (gl, g2)pld(g2) + L 2l2 (?l ) q2) P 2d{ q 2) 

Here components of L-tensor L* (q) are the unknown functions on the configuration space Q. 
This L-tensor gives rise to a torsionless tensor field N of type (1,1) on T*Q, which acts on the 
1-form of fibcred coordinates (p, q) as 

> Nq:={seq(d(q[k])=add(L[k,j] (q)*d(q[j]) , j=l . .n) ,k=l . .n)} : 

> Np:={seq(d(p[k])=add(L[j,k] (q)*d(p[j]) 

> - p| I j*add((diff (L[j ,i] (q),q[k]) -diff (L[j,k] (q) ,q[i] ))*d(q| |i),i=l. .n), 

> j=l. .n) ,k=l. .n)>: 

At the second step we introduce canonical Poisson tensor & 

> ed:=array(identity, l..n,l..n): u:=array( sparse, 1 . .n, 1 . .n) : 

> P : =linalg [stackmatrix] (linalg [augment] (u,ed) , 

linalg [augment] (-ed,u)) ; 

10" 

1 

-1 
0-100 

and canonical Poisson brackets 



P 
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> PB : =proc(f ,g) options operator, arrow: 

> add( diff (f ,p| |i)*diff (g,q| |i)- diff (f ,q| |i)*diff (g,p| I i) , i=l..n) 

> end: 

It allows us to calculate the vector fields Xt — &dT(p, q) and Xy = 0^dV{q). They are equal 
to 

> dT:=array(l . .2*n) : 

> for i from 1 to 2*n do dT [i] : =dif f (T(var) , var [i] ) : end do: 

> X_T:=evalm(P&*dT) ; 



XT:= | -^jT(ql, q2,pl,p2), -^T(ql, q2,pl,p2), 

- (^7 T (^> q2,pl,p2)), -(^ T (^: q2,pl,p2)) 

and 

> dV:=array(l . .2*n) : 

> for i from 1 to 2*n do dV [i] : =dif f (V(q) , var [i] ) : end do: 

> X_V:=evalm(P&*dV) ; 



XV 



The equations with the vector field Xy will be simpler due to its zero components than ones 
with Xt and, therefore, we start with the second equation (|2.12l) in order to explain some 
features of the Maple procedures. Namely, substituting the vector field Xy and 1-form 9 into 
the ±x v d9, one gets 



> idT:=wcollect(hook(d(theta) ,X_V)) : 



icff:=( - L 2A (ql,q2)(—V(ql,q2))(ql,q2,pl,p2) 

- L ltl (ql, rz2)(A V (gl, q 2))(ql, q2, pi, p2))d(?l) 
+ ( - L 2>2 (ql, ,2)(A V ( g l, q2))(ql, q2, pi, p2) 

- L 1>2 (ql, ,2)(A V ( g l, q2))(ql, q2, pi, p2))d(g2) 

It is easy to see that this expression contains derivatives of the potential V(q), which formally 
depend on all the variables (q,p). In fact these derivatives depend on the variables q only. It is 
an unpleasant feature of the Maple procedure hook from the package liesymm. 

In order to get the correct expression for ±x v d9 we have to use the special substitution 

> Trans :={seq(diff (V(q),q| I i) (var)=dif f (V(q) ,q| | i) ,1=1. .n)}: 

> idT : =subs (Trans , idT) ; 



idT:=( L 2tl (ql, q2) — V(ql, q2) - L 1A {ql, q2) — V(ql, g2))d(gl) 

+ ( - L 2 , 2 (ql, q2) -^V(ql, q2) - L 1>2 {ql, q2) A V(gl, g2))d( 5 2) 

The equation l|2.12|l are satisfied identically and coefficients for independent 2-forms in l|2.12|l 
vanish. It gives rise to the following system of equations 
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> SysEq: =annul(d(idT-V(q) *d(sigma) ) , [var] ): nops(SysEq); 

At ft = 2 equation i|2.12[) generates only one partial differential equation, while at n = 3 the 
system SysEq consists of three equations and so on. 
We keep these equations in a special list 

> ListEq:=NULL: 

> for i from 1 to nops (SysEq) do 

> ListEq:=ListEq,lhs(SysEq[i] ) ; 

> end do : 

Now we pass on the remaining equation p.lljl . As above we have to simplify this equation 
through the additional substitution for derivatives after the application of the hook operator 

> Eq: =wcollect (d(hook(d(theta) ,X_T) -T(var) *d(sigma) ) ) : 

> Trans :={seq( dif f (T(var) ,var [i] ) (var)=diff (T(var) ,var [i] ) ,i=l . .2*n)}: 

Expanding equation (|2.11|l in the basis of 2-forms one gets one more system of equations 

> SysEq: =anmil (subs (Trans ,Eq) , [var] ) : nops(SysEq); 

At n = 2, 3 equation l|2.12|l gives rise to six and fifteen equations respectively. 

By adding these equations to earlier prepared list of equations ListEq one gets a complete 
set of algebraic and partial differential equations on the components of L-tensor. 

> for i from 1 to nops(SysEq) do ListEq:=ListEq,lhs(SysEq[i] ) ; end do: 

> NumEq : =nops ( [ListEq] ) ; 

NumEq := 7 

At n — 2, 3 one gets seven and eighteen equations on four and nine functions Lj(q) respectively. 

At the third step we have to define the Hamilton function. As an example we consider the 
Henon-Heiles integrable model with the Hamiltonian 

h = P \ + v\ + + ql) + qfa + >l ■ (3.13) 

In this case kinetic energy T and potential V{q) are given by 

> Tn:=add(p| |i~2,i=l. .n) ; 

> Vn:=l/2*(ql~2+q2~2)+ql"2*q2+2*q2"3; 

Substituting these parts of the Hamilton function into the prepared above list of equations 
ListEq one gets a set of polynomial equations of the second degree in the momenta p, which 
must be identically satisfied for all admissible values of variables p. It means that coefficients 
for the second, first and zeroth power of pi vanish. All these coefficients form new system of 
algebraic and partial differential equations on the functions (q) . 

> System := NULL: 

> for i from 1 to NumEq do 

> u:=expand(dvalue(subs({ T(var)=Tn, V(z)=Vn }, ListEq [i] ))) : 

> System:=System,coef f s(u,{p}) : 

> end do : 

> NumSys : =nops ({System}) ; 

NumSys := 13 

At n — 2, 3, 4 the resulting system of equation consists of 13, 51 and 136 equations. Even though 
some of these equations may be dependent, it is convenient to use them all simultaneously. Of 
course, for a generic potential this system is heavily overdetermined and has only the trivial 
solution, which means that the Hamilton function is non-separable. 

In the case of polynomial potentials we can easy solve this system by means of the standard 
Maple procedure pdsolve 



G 



> Ans :=pdsolve( {System}, {seq(seq(L[i, j] (z) , i=l . .n) , j=l . .n)} ) 

which allows us to get the following answer 



Ans := {L hl (ql, q2) 



3Ci 



+ C 2 , L 1>2 (ql, q2) 



L 2t2 (ql, q2)=C x q2 + C 2 , L 2 ,i{ql, q2) 

depending on two arbitrary constants Ci,2- 

Substituting this answer into the L-tensor one gets 

> Trans :={seq(seq(L[i,j] (q)=L[i,j] ,i=l. .n) ,j=l. .n)}: 

> Ans : =simplify(map2 (subs , Trans , Ans) ) : 

> L:=array(l . .n, 1 . .n) : 

> L :=map2( subs, Ans, evalm (L) ) ; 



L : = 



3 r 



Co 



c iq i 



c iq i 



CI q2 + C* 2 



The eigenvalues of this L-tensor are the separation coordinates 
> Q : =linalg [eigenvalues] (Ln) ; 



Ciqi 
2 

Ci ql 



} 



Q-= 3 -^ + c 2 + % 



9 C\ - 24 C\ q2 + 16 C? q2 2 + 16 C? ql '' 



2 „1 2 



3d ri C lQ 2 yj9Cl-24Cfq2 + WCfq2' + WClql 
— h C 2 H 



Thus one gets translated parabolic coordinates, which are the separation coordinates for the 
Henon-Heiles model [7]. It takes less then a minute of the computer time. 

More explicitly, we obtain a family of equivalent separable coordinate systems labeled by 
Ci, 2 - Recall, that two separable systems are called equivalent if the corresponding separated 
solutions of the Hamilton- Jacobi equation generate the same Lagrangian foliation of T* Q. 

Now let us consider construction of the corresponding integrals of motion in the framework 
of the Riemannian geometry. Following Benenti we have to construct polynomials a m 

> for i from to n do 

> sigma[i] : =coef f (linalg[det] (lambda*ed-L) , lambda, n-i) : 

> end do : 

and basis of Killing tensors K m (I2.4I) 

> for m from 1 to n-1 do 

> K| |m:=evalm(add( sigma[m-k] *L~k,k=0 . .m) ) : 

> end do : 

Remind, that in our case g y = 1 and, therefore, L y = L*. 

In addition to K m equations (|2.5() consist of exterior derivatives on the unknown potentials 
dV m , which are defined by 

> for m from to n-1 do 

> dV I |m: =array (1 . .n) : for i from 1 to n do 

> dVl |m[i] :=diff(V| |m(q),q[i]): 

> end do : 

> end do : 

> dV0:=map2(subs,V0(q)=Vn,dV0) : 
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Now we can solve equations l|2.5|l and determine integrals of motion H m Q2.6JI 

for m from 1 to n-1 do 

> eqV:=evalm(dV| |m-K| |m&*dV0) : 

> Ans:=pdsolve({seq(eqV[i] , i=l . .n) } , V I I m(q) ) : 

> H[m] :=add(add( K| |m[i,j]*p| |i*p| I j ,i=l. .n) ,j=l. .n)+subs(Ans,V| |m(q)) : 

> end do : 

We can show the Hamilton function 

> H[0] :=Tn+Vn; 

H := pi 2 + p2 2 + ^— + ^j- + ql 2 q2 + 2 q2 3 
and prove that our integrals of motion H m commute with the Hamiltonian Hq 

> for m from 1 to n-1 do 

> ZERO:=simplify(PB(H[0] ,H[m] )) ; 

> end do ; 

ZERO := 



Remark: At n < 10 and for the polynomial potentials we can solve overdetermined system of 
equation System with the standard procedure pdsolve on a standard personal computer in the 
reasonable time. At n > 10 we have to use special computers or special symbolic software for 
solving heavily overdetermined systems of equations. 



4 Examples 

In this section we present some examples which illustrate various aspects of the proposed code. 



4.1 The anharmonic oscillator 

Configuration space Q is the n-dimensional Euclidean space R ra with cartesian coordinates q 
such that metric takes the form ds 2 = ^ gijdqidqj — ^2 dqi 2 . Let us consider the anharmonic 
oscillator with the Hamiltonian 

n n / n \ 2 

i=i i=i V i=i / 

We have to substitute this Hamiltonian instead of (|3.13|l using commands 

> Tn:=add(p| |i~2,i=l. .n) ; 

> Vn:=add(a| |i*q| |i~2,i=l. .n)+(add(q| |i~2,i=l. .n))~2; 

In is easy to calculate solution of the complete system of equations System obtained from 
equations (|2.11|) and I2.12|l 

L l 3 (q) = (Ci + di - a n )Sij + C 2 qtqj 

Here the solution found by Maple is presented without any further simplification. 

At C\ = a n and C2 = 1 the eigenvalues Qi of the corresponding tensor L are defined by 
the equation 

det ( £ - A ) = _i + y^___-n A -^ (4i4) 
nr =1 (A-«o- + ^A- ai iiA-ar (4 - i4j 

This is the well-known relation which defines standard elliptic coordinates in K n . 
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These coordinates were introduced by Jacobi in a note in Crelle's Journal ^I] • A thorough 
discussion of its general properties as well as of its use for separation of variables in the Hamilton- 
Jacobi equation can be found in his lecture notes [Q. 

Now let us consider construction of integrals of motion in the involution in framework of 
the bi-Hamiltonian geometry. At first we determine action of the recursion operator N on the 
1-forms dp and dq using given tensor L 

> Nqn:=wcollect(simplify(map2(subs,Ans,Nq))) : 

> Npn:=wcollect(simplify(map2(subs,Ans,Np))) : 

Then we calculate symmetric polynomials a m on eigenvalues of L 

> for m from to n do 

> sigma[m] : =coef f (linalg [det] (Ln-lambda*ed) , lambda, n-m) : 

> end do : 

At the final step we construct and solve recursion relations 12.7( 1 

> unassign(H): H[n+1]:=0: H[0]:=Tn+Vn; 

> for m from n-1 by -1 to 1 do 

> Eq:=wcollect(dvalue( d(H[m+l] )- 

> subs(Nqn union Npn,d(H [m] (var) ) ) -sigma[m+l] *d(H [0] ) ) ) : 

> Sys : =annul (Eq, [var] ) : 

> Ansl : =pdsolve (Sys ,H [m] (var) ) ; 

> H[m] :=subs (Ansl, H[m] (var)) : 

> end do : 

In addition we can check the Poisson brackets relations between integrals of motion 

> for m from 1 to n-1 do 

> ZERO:=simplify(PB(H[0] ,H[m] )) ; 

> end do ; 



4.2 The Euler two centers problem 

Let us consider the Euler problem of planar motion in a force field of two attracting centers 
with the Hamiltonian 



H =y{ 



Pi 



a i 



o 2 



y/(qi - c) 2 + q\ y / (q 1 + c) 2 + q\ 



In contrast with the previous examples potential V{q) is an algebraic function instead of a 
polynomial one. 

Solution of the overdetermined system of equations obtained from 12.1111 and 1(2.12(1 is equal 

to 



L = 



(q{ - c 2 )C 2 + d q x q 2 C 2 



<7l<?2C2 



q 2 2 C 2 + d 



At Ci = 1/2 c 2 and C 2 = \ the eigenvalues of the matrix L are the standard elliptic coordinates 



_ (A-Qi)(A-Q 2 ) 

A + c 2 /2 A-c 2 /2 (A-c 2 /2)(A + c 2 /2) 



qi 



ql 



4.3 The Toda lattice 

As above, configuration space Q is the n-dimensional Euclidean space IR ra with cartesian coor- 
dinates q. The Hamiltonian of the periodic Toda lattice is given by 

_^ n n 

H = - ^Pi 2 +^2exp(q, - q i+ i) , q n+1 = q x . (4.15) 
i=i i=i 
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At 7i = 2 there exists the unique solution of the complete system of equations System 



L 



C\ C2 
C2 C\ 



Of course, eigenvalues of this matrix are not the separation variables. However, in this case the 
number of free parameters is equal to the dimension of the Riemannian manifold Q. It allows 
us to construct the separation coordinates too (see [7j). Namely, let us diagonalize matrix L 



L = l/- 1 diag(Ci + C 2 , Ci - d)V, 



V = 



1 1 

-1 1 



It is easy to see that the separation coordinates Q are the following cartesian coordinates 



Q = Vq, 



h = Ql + Q2, Q2 = qi + Q2 ■ 



At n = 3 solution of the equations System is the trivial o 



L : 



C\ C2 C2 
C 2 C\ C2 
C2 C2 C\ 



The number of free parameters is less than the dimension of the Riemannian manifold. It means 
that this Hamiltonian does not separable through point transformations. However, we have to 
underline that this Hamiltonian admits separation of variables in a wider class of canonical 
transformations of the whole phase space. 



4.4 The Neumann system 

The unit sphere S™ is the Riemannian subspace of K n+1 whose points have cartesian coordinates 
x = (xi, . . .,x n +i) satisfying \x\ = \/ (x, x) = 1. 

The Neumann system is a well-known and very much studied mechanical system, which 
describes particle moving on the sphere under the influence of a quadratic potential V = 
h ^2 dixf. Due to the form of the constraint and of the potential, it is quite natural to perform 
all the computations in the local coordinates = x\, i = 1, . . . , n. In these coordinates the 
Hamiltonian of the Neumann system is given by H = T + V where 



i=l 



-22qiqjPiPj, 

i<j 



V 



1 

- ^(oj - a n+1 )qi 



Here pi are momenta conjugated to qi and a\, . . . , a ra +i are arbitrary parameters |l()j . We can 
insert this Hamiltonian in the Maple code with the commands 

> Tn:=2*add(q| |i*(l-q| |i)*p| |i~2,i=l. .n) 

> -4*add(add(q| |i*q| I j*pl |i*pl I j ,j=i+l. .n) ,i=l. .n) ; 

> Vn:=l/2*add((a| |i-a| I (N+l))*q| 1 1 , 1=1 . .n) ; 

At n = 2 the unique solution of the complete system of equations System looks like 

qi(a 3 - 01) + ai - a 2 



L 



O2- 



a 2 - a 3 
92 (gi - 03) 
a 2 - a 3 



■ C x qi Ci 

q 2 Ci + C 2 



Recall, that L^(q) are components of the tensor L, which is symmetric with respect to the 
metric G. For the non-flat metrics the L-matrices obtained in Maple are non-symmetric. 
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If C\ = a 3 — a 2 and C 2 = a 2 the eigenvalues Q\. 2 of this L-tensor are zeroes of the following 
function 

det(L-A) qi q 2 l-<?i-<?2 

e(A) = n n 77 7 = T + T + 7 

In the redundant coordinates Xi one gets the following definition of the separation coordinates 



— — . with > Xj = 1 . 



These are well-known defining relations for the spheroconical or elliptic-spherical coordinates. 
At n > 2 one gets the same relations too. 

In contrast with the previous examples in this case wehave non- trivial metric and (q) 7^ 
L 13 {q) . In order to get integrals of motion it is more convenient to use recursion relations l|2.7|l . 
The corresponding part of code looks like 

> Nqn:=wcollect(simplify(map2(subs,Ans,Nq))) : 

> Npn:=wcollect(simplify(map2(subs,Ans,Np))) : 

> for m from to n do 

> sigma [m] : =coef f (linalg [det] (Ln-lambda*ed) , lambda, n-m) : 

> end do : 

> unassign(H): H[n+1]:=0: H[0]:=Tn+Vn; 

> for m from n-1 by -1 to 1 do 

> Eq:=wcollect(dvalue( d(H[m+l] )- 

> subs(Nqn union Npn,d(H [m] (var) ) )-sigma[m+l] *d(H[0] ) ) ) : 

> Sys : =annul (Eq, [var] ) : 

> Ansl : =pdsolve (Sys ,H [m] (var) ) ; 

> H[m] :=subs (Ansl, H[m] (var)) : 

> end do : 

Of course, we can check the Poisson brackets relations between these integrals of motion 

> for m from 1 to n-1 do 

> ZERO:=simplify(PB(H[0] ,H[m])) ; 

> end do ; 

At 71 = 2, 3 all the calculations take about one minute and eight minutes respectively. 

Remark: In contrast with [TO] we directly solve equations (|2.8I2.10|) without any additional 
assumptions about the affine structure of the solutions. It allows us to prove that there is one 
unique solution only. 

4.5 The Jacobi-Calogero inverse-square model 

Configuration space Q is 3-dimensional Euclidean space R 3 with cartesian coordinates q = 
(<7i, qi , (73). The Hamilton function is given by 

3 

H = ^2 Pi + {qi - q-z)' 2 + (92 - 93)^ + (93 - 9i)~ 2 • 

i=l 

For this Hamiltonian solution of the complete system of equations System depends on four 
arbitrary parameters 



L = 



q\C x + 2 qi C 2 + C 3 qxqiCx + (91 + q 2 )C 2 + C 4 9193C1 + (gi + q 3 )C 2 + C 4 

9l92Ci + (ft + q 2 )C 2 + C 4 q\C x + 2q 2 C 2 + C 3 g 2 q 3 C*i + (q 2 + q 3 )C 2 + C 4 

9193C1 + (91 + q 3 )C 2 + Ci q 2 q 3 Ci + (q 2 + q 3 )C 2 + C" 4 qjCi + 2q 3 C 2 + C 3 
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The eigenvalues of this L-tensor are the separation coordinates 



> Q : =linalg [eigenvalues] (L) ; 

The number of free parameters Ck is more than the dimension of the Riemannian manifold. 
It means that this integrable system is degenerate or superintegrable system, which admits 
separation of variables in some different coordinate systems. Namely, at the different values of 
the parameters C\, . . . , C4 eigenvalues of the tensor L are oblate spheroidal, prolate spheroidal, 
spherical, rotational parabolic and circular cylindrical coordinates. The detailed discussion may 
be found in 0G2]. 

4.6 An elliptic egg 

Let us consider the Hamiltonian with the rational potential 



In contrast with the previous examples potential V(q) is a rational function rather than a 
polynomial one. This potential has a singular surface, which is an ellipsoid 



For c > 0, this singular surface is repelling inwards, so that all trajectories starting inside the 
ellipsoid remain there forever. 

At n > 2 the Maple procedure pdsolve can not solve the system of equations System 
obtained from the both equations (|2.11|) and (|2.12|l in the reasonable time. 

In this case we have to solve equation for the geodesic motion (|2.11|) and then the equation 
on potential (|2.12l) . At n < 10 it takes just a few minutes. 

In our case solution of the system of equations generated by the geodesic equation l|2.11[l 
is equal to 



Here coefficients a,Aij,Bij and Cy are arbitrary constants, Aij = Bji and dj — C/j. 

Substituting this solution into the system of equations obtained from H2.12|) one gets a 
system of algebraic equations on the coefficients Ay , Bij and Cy . The standard Maple program 
solve easy solves this system of algebraic equations in few seconds. The eigenvalues of the 
corresponding tensor L are the standard elliptic coordinates (|4.14l) . 

5 Conclusion 

We present the first part of the symbolic software which builds the separation coordinates in 
the Hamilton- Jacobi equation for the L-systems. The second part of this software is devoted to 
building and solving the corresponding separated equations. We leave this part to subsequent 
publications. 

Recall, that Maple file with the first part of the code may be found in the following URL 
http : //www . maplesof t . com/applications/app_center_view . aspx?AID=1686 

The latest version of the Maple file may be obtained from authors. 





L){q) = aq l q j + A^qt + BijQj + C tj . 
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